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ABSTRACT 

A  time-dependent  generalization  of  a  Fourier-ray  method  is  presented  and  tested  for  fast  numerical 
computation  of  high-resolution  nonhydrostatic  mountain-wave  fields.  The  method  is  used  to  model  moun¬ 
tain  waves  from  Jan  Mayen  on  25  January  2000,  a  period  when  wavelike  cloud  banding  was  observed  long 
distances  downstream  of  the  island  by  the  Advanced  Very  High  Resolution  Radiometer  Version  3 
(AVHRR-3).  Surface  weather  patterns  show  intensifying  surface  geostrophic  winds  over  the  island  at  1200 
UTC  caused  by  rapid  eastward  passage  of  a  compact  low  pressure  system.  The  1200  UTC  wind  profiles  over 
the  island  increase  with  height  to  a  jet  maximum  of  —60-70  m  s“^,  yielding  Scorer  parameters  that  indicate 
vertical  trapping  of  any  short  wavelength  mountain  waves.  Separate  Fourier-ray  solutions  were  computed 
using  high-resolution  Jan  Mayen  orography  and  1200  UTC  vertical  profiles  of  winds  and  temperatures  over 
the  island  from  a  radiosonde  sounding  and  an  analysis  system.  The  radiosonde-based  simulations  produce 
a  purely  diverging  trapped  wave  solution  that  reproduces  the  salient  features  in  the  AVHRR-3  imagery. 
Differences  in  simulated  wave  patterns  governed  by  the  radiosonde  and  analysis  profiles  are  explained  in 
terms  of  resonant  modes  and  are  corroborated  by  spatial  ray-group  trajectories  computed  for  wavenumbers 
along  the  resonant  mode  curves.  Output  from  a  nonlinear  Lipps-Hemler  orographic  flow  model  also 
compares  well  with  the  Fourier-ray  solution  horizontally.  Differences  in  vertical  cross  sections  are  ascribed 
to  the  Fourier-ray  model’s  current  omission  of  tunneling  of  trapped  wave  energy  through  evanescent  layers. 


1.  Introduction 

When  suitable  environmental  conditions  exist,  flow 
over  mountains  generates  quasi-stationary  gravity 
waves  that  can  propagate  obliquely  away  from  the  par¬ 
ent  orography.  These  waves  can  have  important  effects 
on  other  atmospheric  processes.  Air  parcels  advected 
through  these  waves  experience  rapidly  oscillating  adia¬ 
batic  cooling  and  heating  that  can  have  strong  net  in¬ 
fluences  on  cloud  physics  and  associated  chemistry  in 
both  the  troposphere  and  stratosphere  (e.g.,  Jensen  et 
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al.  1998;  Fueglistaler  et  al.  2003).  These  influences 
sometimes  take  the  form  of  spectacular  wave-banded 
cloud  displays  that  are  visible  from  space  (e.g.,  Fritz 
1965;  Gjevik  and  Marthinsen  1978;  Burroughs  and  Lar¬ 
son  1979;  Sharman  and  Wurtele  1983;  Mitchell  et  al. 
1990;  Worthington  2001).  Mountain  waves  grow  in  am¬ 
plitude  with  altitude  and  can  break,  generating  drag 
forces  that  affect  the  synoptic-scale  circulation  and  tur¬ 
bulence  that  mixes  chemical  species  (e.g.,  Kim  et  al. 
2003). 

In  addition  to  these  atmospheric  effects,  mountain 
waves  present  hazards  to  aviation.  Severe  structural 
damage  and  injuries  to  passengers  and  crew  can  occur 
when  aircraft  fly  through  severe  clear- air  turbulence 
(CAT)  produced  by  mountain-wave  breaking  (e.g., 
Bacmeister  et  al.  1994;  Ralph  et  al.  1997;  de  Villiers  and 
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van  Heerden  2001).  Furthermore,  sudden  changes  in 
flight  altitude  caused  by  mountain  waves,  either 
through  wave-induced  CAT  and/or  short-wavelength 
vertical  wave  displacements,  are  also  an  important  issue 
given  recent  enaction  of  reduced  vertical  separation 
minima  (RVSM)  for  commercial  air  traffic.  The  large- 
amplitude  mountain-wave  event  over  Colorado  docu¬ 
mented  by  Lilly  (1978)  provides  a  vivid  illustration  of 
both  of  these  hazards,  causing  a  commerical  airliner  to 
drop  —4000  ft  in  a  little  over  1  min  while  simulta¬ 
neously  undergoing  severe  airframe  buffeting  due  to 
wave-induced  CAT. 

Thus,  mountain  wave  events  are  important  to  fore¬ 
cast,  yet  the  spatial  scales  of  these  waves  and  the  sub¬ 
wavelength  instabilities  that  lead  to  breaking,  drag,  and 
turbulence  generation  are  generally  too  short  for  op¬ 
erational  numerical  weather  prediction  (NWP)  models 
to  resolve  fully  (e.g.,  Benoit  et  al.  2002;  Smith  2004). 
For  example,  subgrid-scale  gravity  wave  drag  must  be 
parameterized  in  NWP  and  climate  models  (e.g.,  Kim 
et  al.  2003).  Until  operational  NWP  model  resolutions 
improve,  alternative  forecasting  algorithms  for  moun¬ 
tain  waves  must  be  developed. 

Ongoing  efforts  in  this  area  have  been  undertaken  at 
the  Naval  Research  Laboratory  for  over  a  decade  in 
developing  the  MountainWave  Forecast  Model 
(MWFM;  Eckermann  et  al.  2004,  2006).  Instead  of 
simulating  the  fully  nonlinear  discretized  Navier- 
Stokes  equations,  as  in  a  traditional  NWP  model,  the 
MWFM  approach  uses  ray  methods  to  simulate  the 
generation,  propagation  and  breakdown  of  mountain 
waves  within  the  large-scale  environment  specified  by 
operational  NWP  model  output.  Since  ray  solutions  are 
easy  to  interpret,  computationally  fast,  and  valid  within 
a  broad  variety  of  atmospheric  environments,  they  are 
attractive  as  potential  forecasting  algorithms. 

The  first  version  of  the  MWFM  (MWFM-1)  tested 
these  ideas  using  a  hydrostatic  two-dimensional  spatial- 
ray  formulation  (Bacmeister  et  al.  1994).  The  MWFM-2 
was  extended  to  use  three-dimensional  spatial-ray 
equations  governed  by  a  nonhydrostatic  dispersion  re¬ 
lation  including  rotation  (Eckermann  and  Preusse 
1999).  Both  models  employ  an  idealized  ridge  decom¬ 
position  of  the  earth’s  topography  to  define  the  major 
terrain  features  relevant  for  generation  of  waves  near 
the  surface.  MWFM  forecasts  have  been  used  for  over 
a  decade  now  to  direct  National  Aeronautics  and  Space 
Administration  (NASA)  research  aircraft  away  from 
hazardous  mountain-wave-induced  turbulence  and  into 
regions  where  nonbreaking  waves  generate  cloud,  tasks 
for  which  it  has  repeatedly  shown  skill  (e.g.,  Bacmeister 
et  al.  1994;  Eckermann  et  al.  2004,  2006).  MWFM-2 
hindcasts  of  stratospheric  wave  amplitudes  have  been 


validated  globally  using  new  data  from  stratospheric 
research  satellites  (Eckermann  and  Preusse  1999;  Jiang 
et  al.  2002,  2004),  as  well  as  regionally  using  various 
suborbital  measurements  (e.g.,  Hertzog  et  al.  2002; 
Eckermann  et  al.  2006).  MWFM  forecasts  of  strato¬ 
spheric  mountain-wave  turbulence  were  utilized  exten¬ 
sively  by  the  U.S.  Air  Force  during  Operations  Endur¬ 
ing  Freedom  and  Iraqi  Freedom  (Eckermann  2002), 
and  since  2004  MWFM-2  has  been  run  operationally  at 
the  Air  Force  Weather  Agency.  MWFM-2  hindcasts 
have  also  been  used  in  a  variety  of  research  applica¬ 
tions,  such  as  the  role  of  mountain  waves  in  accelerating 
ozone  loss  chemistry  in  the  Arctic  winter  stratosphere 
(e.g.,  Carslaw  et  al.  1999;  Pierce  et  al.  2003;  Svendsen  et 
al.  2005;  Mann  et  al.  2005). 

Thus,  ray  methods  are  now  an  accepted  approach  to 
forecasting  mountain  waves  (Eckermann  et  al.  2004). 
Yet  the  current  MWFM  ray  algorithms  still  contain  a 
number  of  significant  simplifications  and  shortcomings. 
First,  they  use  a  one-dimensional  height-dependent 
form  for  the  wave  action  equation  that  does  not  include 
the  effects  of  horizontal  geometrical  spreading  of  the 
rays  on  the  wave  amplitude  evolution  (e.g.,  Shutts 
1998).  Second,  they  use  a  spatial  formulation  for  the  ray 
solutions,  which  leads  to  caustic  singularities  that  are 
not  practical  to  correct  (e.g.,  Broutman  et  al.  2001, 
2002,  2004).  Third,  they  use  idealized  ridge  databases 
for  their  source  functions  that  do  not  capture  the  full 
spectrum  of  waves  radiated  by  flow  over  realistic  to¬ 
pography. 

To  fully  explore  and  exploit  the  capabilities  of  ray 
methods  for  mountain-wave  forecasting,  over  the  past 
several  years  we  have  progressively  developed  im¬ 
proved  ray  algorithms  that  reduce  or  eliminate  these 
and  other  weaknesses  in  the  current  MWFM  ray  equa¬ 
tions  (Broutman  et  al.  2001,  2002,  2003,  2004,  2006). 
This  has  led  to  a  new  ray  algorithm  that  we  refer  to  here 
as  the  Fourier-ray  method,  because  it  involves  Fourier- 
synthesized  (rather  than  spatially  synthesized)  ray  so¬ 
lutions.  The  method  has  recently  been  reviewed  inter 
alia  by  Broutman  et  al.  (2004)  and  the  version  of  it  that 
we  use  here  is  described  in  section  2a. 

This  method  shows  promise  as  a  potential  next- 
generation  dynamical  ray  core  for  the  MWFM,  since 
the  idealized  solutions  derived  to  date  alleviate  or 
eliminate  all  of  the  aforementioned  weaknesses  of  the 
corresponding  spatial  ray  solutions.  Specifically,  the  al¬ 
gorithm  incorporates  arbritrary  topographic  forcing  at 
the  lower  boundary,  models  both  trapped  and  free- 
propagating  waves,  includes  horizontal  geometric 
spreading  in  its  wave  action  solutions,  and  systemati¬ 
cally  corrects  the  caustic  singularities  found  in  spatial- 
ray  solutions.  To  date,  however,  we  have  applied  these 
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Fig.  1.  Three-dimensional  shaded  surface  plot  of  topographic  elevation  for  the  island 
of  Jan  Mayen  from  the  GTOPO30  digital  elevation  database.  All  three  length  dimen¬ 
sions  (in  km)  as  well  as  longitudes  and  latitudes  are  displayed. 


new  algorithms  only  to  idealized  mountain-wave  prob¬ 
lems  that  use  analytical  functions  for  the  obstacle  shape 
and  background  wind  profile. 

Here  we  take  our  next  step  in  assessing  the  potential 
of  the  Fourier-ray  method  for  forecasting  by  extending 
and  applying  these  algorithms  to  a  more  realistic  case 
study.  The  major  extensions  here  are  incorporation  of 
realistic  topography  and  real  background  wind  and 
temperature  profiles  from  both  radiosonde  observa¬ 
tions  and  a  meteorological  analysis  issued  by  a  data 
assimilation  system.  We  concentrate  here  on  nondissi- 
pative  wave  solutions,  deferring  developments  in  esti¬ 
mating  locations  of  wave-induced  CAT  to  later  studies. 

Our  modeling  focuses  on  Jan  Mayen  (71°N,  8.4°W), 
a  small  island  in  the  North  Atlantic  whose  topography 
is  dominated  to  the  north  by  Mount  Beerenberg,  a 
quasi-circular  volcanic  mountain  of  width  —5-10  km 
and  a  peak  elevation  of  —2270  m  (see  Fig.  1).  We  seek 
to  “hindcast”  wavelike  patterns  observed  over  the  is¬ 
land  in  satellite  cloud  imagery  on  25  January  2000.  Pre¬ 
vious  analysis  and  linear  modeling  of  some  earlier  sat¬ 
ellite  cloud  images  over  Jan  Mayen  (Gjevik  and  Mar- 
thinsen  1978;  Simard  and  Peltier  1982)  associated  cloud 
banding  here  with  three-dimensional  trapped  lee  waves 


forced  by  flow  over  Mount  Beerenberg  that  radiated 
into  much  larger  atmospheric  volumes  downstream. 
Such  waves  present  an  excellent  test  case  for  our  new 
algorithm’s  vertical  reflection  and  high-resolution  fore¬ 
casting  capabilities,  since  Jan  Mayen’s  topography 
would  not  be  adequately  resolved  by  current  opera¬ 
tional  NWP  systems.  Our  results  here  are  used  to 
benchmark  the  initial  performance  of  this  new  ray  code 
and  to  target  areas  for  further  development  for  future 
forecasting  applications. 

2.  Models 

a.  The  Fourier-ray  model 

Our  primary  tool  is  a  flexible  numerical  implemen¬ 
tation  of  a  Fourier-ray  method.  The  name  “Fourier 
ray”  refers  to  Fourier-synthesized  ray  solutions:  that  is, 
the  ray  solutions  are  computed  in  a  Fourier  domain  and 
then  superimposed  by  inverse  Fourier  transform  to  give 
a  spatial  solution.  The  method  is  described  in  Brout- 
man  et  al.  (2002,  2003,  2006).  During  earlier  stages  of  its 
development,  different  names  were  used  for  the  same 
basic  method  (e.g.,  Maslov’s  method,  a  simplified  Fou- 
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rier  method).  The  reasons  for  the  name  changes  are 
discussed  in  Broutman  et  al.  (2006). 

The  ray  solution  in  the  Fourier  domain  is  expressed 
as  a  function  of  k,  I,  z,  t,  where  k  and  I  are  the  horizontal 
wavenumbers,  z  is  height,  and  t  is  time.  The  advantage 
of  solving  in  the  Fourier  domain  is  that  k  and  /  are 
constant  along  the  ray  in  a  horizontally  uniform  back¬ 
ground,  as  assumed  here.  Along  with  a  simple  ray- 
based  treatment  of  the  time  dependence  of  the  solu¬ 
tions  (Broutman  et  al.  2006),  the  problem  is  effectively 
reduced  to  a  one-dimensional  calculation  in  z.  This  is  a 
major  simplification  for  numerical  ray  tracing,  espe¬ 
cially  for  the  wave  amplitude  calculations  and  the  cor¬ 
rection  of  caustics. 

Broutman  et  al.  (2002,  2003)  derived  steady-state 
Fourier-ray  solutions  for  hydrostatic  and  nonhydro¬ 
static  waves,  respectively.  Wave  transience  was  incor¬ 
porated  into  the  Fourier-ray  method  by  Broutman  et  al. 
(2006).  We  use  that  transient  formulation  here  to  aid 
direct  comparisons  with  output  from  nonlinear  numeri¬ 
cal  models  at  finite  times.  Besides  being  helpful  for 
estimating  setup  times  for  the  wave  field,  the  transient 
formulation  has  two  computational  advantages.  First, 
there  are  no  resonant  singularities  in  the  transient  so¬ 
lution  for  the  trapped  waves,  as  there  are  in  the  steady- 
state  solution.  (The  steady-state  singularities  can  be  re¬ 
moved  by  adding  a  small  imaginary  component  to  the 
wavenumber  or  frequency,  but  this  artificially  damps 
wave  amplitudes.  The  Fourier-ray  solutions  presented 
here  are  nondissipative.)  Second,  the  transient  wave 
field  is  more  limited  in  spatial  extent  than  the  steady- 
state  solution,  which  is  the  longtime  limit  of  the  tran¬ 
sient  solution.  The  transient  solution  can  thus  be  rep¬ 
resented  with  fewer  computational  grid  points. 

We  consider  linear  three-dimensional  mountain 
waves  radiated  from  realistic  topography  into  arbitrary 
vertical  profiles  of  winds  and  stability.  The  mountain 
waves  are  stationary,  with  zero  frequency  (co  =  0)  rela¬ 
tive  to  the  ground.  For  a  horizontal  background  wind 
U  =  (U,V,0),  the  intrinsic  frequency  isw=co  —  k*U  = 
—kU  —  IV,  where  k  =  {k,  I,  m)  is  the  wavenumber 
vector.  The  background  wind  and  buoyancy  frequency 
N  are  assumed  to  depend  on  height  but  not  on  hori¬ 
zontal  position  or  time. 

The  ray  tracing  performed  here  is  based  on  a  gravity 
wave  dispersion  relation  of  the  form 

(h  =  k^N/{k\  +  (1) 

where  k^  =  {k^  +  We  ignore  the  effects  of  the 

earth’s  rotation  and  compressibility  in  (1)  for  simplicity: 
they  are  included  in  other  versions  of  the  code. 

We  define  w{k,  I,  z,  t)  to  be  the  ray  solution  in  the 


Fourier  domain  for  the  vertical  velocity  associated  with 
the  mountain  waves.  The  corresponding  spatial  solu¬ 
tion  w{x,  y,  z,  t)  is  obtained  by  inverse  Fourier  trans¬ 
form: 

r  oo  r  oo 

w{x,  y,  z,t)  =  \  I  w{k,  /,  z,  dk  dl. 

J  — oo  J  — oo 

(2) 

The  Fourier-ray  solution  for  w(k,  I,  z,  t)  is  given  in  the 
appendix,  based  on  derivations  in  Broutman  et  al. 
(2006).  We  emphasize  that  w(x,  y,  z,  t)  is  not  the  same 
as  the  spatial-ray  solution  for  the  vertical  velocity, 
which  is  obtained  from  the  stationary  phase  limit  of  the 
inverse  Fourier  transform  in  (2).  The  distinction  is  im¬ 
portant  because  the  stationary  phase  approximation, 
but  not  w,  breaks  down  at  caustics  in  the  spatial  domain 
(see,  e.g.,  Shutts  1998;  Broutman  et  al.  2001). 

This  Fourier-ray  code  ingests  arbitrary  vertical  pro¬ 
files  of  wind  and  temperature,  which  are  interpolated 
onto  a  vertical  grid  sufficiently  fine  to  evaluate  the  fol¬ 
lowing  phase  integrals  accurately: 

z  Czt  Cz 

mdz,  \  mdz,  \m\dz,  (3) 

0  J  z  J  zt 

where  z^  is  the  height  of  the  turning  point  for  a  given  (k, 
/).  The  first  integral  is  used  for  vertically  propagating 
waves,  the  second  integral  is  used  for  vertically  trapped 
waves  at  heights  below  their  turning  point  z^,  and  the 
third  integral  is  used  for  vertically  trapped  waves  at 
heights  above  their  turning  point  (see  the  appendix). 
An  absolute  value  appears  in  the  third  integral  because 
m  is  imaginary  above  z^.  These  integrals  are  computed 
here  numerically  by  the  trapezoidal  rule,  using  a  verti¬ 
cal  grid  size  Az  =  0.25  km.  Comparisons  with  solutions 
using  smaller  Az  showed  that  this  grid  size  was  ad¬ 
equate. 

The  code  uses  a  finite  Fourier-series  approximation 
to  the  inverse  Fourier  transform  (2).  In  the  cases  pre¬ 
sented  here,  512  wavenumber  values  were  chosen  for  k 
and  /,  corresponding  to  a  spatial  grid  of  1024  by  1024  in 
X  and  y,  and  a  horizontal  grid  spacing  of  Ax  =  Ay  = 
1  km. 

b.  Nonlinear  orographic  flow  model 

To  assess  more  rigorously  the  accuracy  of  the  Fou¬ 
rier-ray  solutions  to  follow,  we  also  performed  some 
companion  simulations  using  a  fully  nonlinear  numeri¬ 
cal  model.  The  model  in  question  solves  the  incom¬ 
pressible  nonlinear  Navier-Stokes  equations  over  to¬ 
pography  using  the  Lipps-Hemler  anelastic  approxima¬ 
tion  (Lipps  and  Hemler  1982;  Nance  and  Durran  1994). 
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The  model  uses  terrain-following  coordinates  (Gal- 
Chen  and  Somerville  1975),  a  centered  second-order 
spatial  discretization,  and  a  modified  leapfrog  time¬ 
stepping  scheme  (Rees  1988).  It  was  run  here  on  a  stag¬ 
gered  Arakawa  C  grid  of  256  points  in  each  horizontal 
direction  and  32  points  in  the  vertical,  with  a  grid  spac¬ 
ing  of  1  km  in  all  directions  and  a  time  step  of  4  s.  A 
free-slip  (frictionless)  lower  boundary  condition  was 
used.  To  absorb  waves  at  heights  ^20  km,  we  imposed 
Rayleigh  damping  with  a  rate  coefficient  that  increased 
linearly  from  zero  at  20  km  altitude  to  4  10“^  s“^  at  the 
32-km  upper  boundary.  To  minimize  reflections  from 
the  lateral  boundaries,  a  radiative  scheme  was  used 
(Miranda  and  James  1992)  as  well  as  some  linearly  in¬ 
creasing  viscous  damping  within  a  10-km-wide  region  at 
each  side  boundary.  Away  from  these  upper  and  side 
boundaries,  the  only  diabatic  terms  were  a  Richardson- 
number-dependent  first-order  turbulent  closure 
scheme  (Lilly  1962)  and  a  small  amount  of  fourth-order 
diffusion  to  suppress  grid-scale  noise  (hyperviscosity 
coefficient  of  10^  m"^  s“^). 

c.  Topography 

Topographic  elevations  h{x,  y)  for  Jan  Mayen  were 
obtained  from  the  United  States  Geological  Survey 
(USGS)  GTOPO30  database.  This  global  dataset  of 
digital  terrain  elevation  has  30  arc  s  resolution,  corre¬ 
sponding  roughly  to  1-km  spatial  resolution,  which  we 
interpolated  linearly  to  our  horizontal  grid  spacing  of 
1  X  1  km^  using  Delaunay  triangulation  and  full  spheri¬ 
cal  to  Cartesian  transformations.  This  model  topogra¬ 
phy  for  Jan  Mayen  is  plotted  in  Fig.  2a  and  can  be  used 
in  the  1  km  X  1  km  Fourier-ray  model  runs. 

For  the  1  km  X  1  km  Lipps— Hemler  model  runs, 
however,  this  1  km  X  1  km  topography  must  be 
smoothed  to  reduce  the  potential  for  unphysical  four- 
point  gridpoint  noise  in  the  simulations  (see,  e.g., 
Davies  and  Brown  2001).  We  achieve  this  using  a  two- 
dimensional  five-point  running  average,  and  the  results 
are  shown  in  Fig.  2b.  In  the  Fourier-ray  model  runs 
shown  here  we  also  use  this  smoothed  topography  to 
facilitate  more  direct  comparisons  with  the  Lipps- 
Hemler  model  output. 

3.  Observational  data 

a.  Advanced  Very  High  Resolution  Radiometer 
Version  3 

The  Advanced  Very  High  Resolution  Radiometer 
Version  3  (AVHRR-3)  is  a  cross-scanning  passive  ra¬ 
diometer,  first  deployed  on  the  National  Oceanic  and 
Atmospheric  Administration  {NOAA-15)  polar-orbit- 


Fig.  2.  Topographic  elevation  contours  h(x,  y)  for  Jan  Mayen, 
(a)  The  unsmoothed  topography  interpolated  to  1  km  X  1  km 
resolution  (as  in  Fig.  1).  (b)  The  topography  after  five-point 
smoothing.  The  contour  interval  is  100  m,  and  the  maximum 
height  of  the  smaller  peak  to  the  lower  left  in  the  plots  is  567  m 
(unsmoothed)  and  359  m  (smoothed). 

ing  satellite.  It  has  six  channels:  three  thermal  infrared 
(IR)  channels  and  three  other  channels  in  the  visible 
and  near  IR.  The  system  acquires  radiances  with  a 
20.32-cm-diameter  telescope  that  scans  cross  track  us¬ 
ing  a  continuously  rotating  mirror.  The  telescope’s  in¬ 
stantaneous  field  of  view  (IFOV)  is  0.0745°,  corre¬ 
sponding  to  surface  footprint  diameters  across  track 
and  along  track  of  1.1  km  X  1.1  km  at  nadir  and  6.2 
km  X  2.3  km  at  the  far  off-nadir  positions.  For  each 
scan,  2048  samples  are  acquired  at  off-nadir  angles  be¬ 
tween  ±55.37°  cross  track  of  the  subsatellite  point,  cor¬ 
responding  to  a  total  horizontal  swath  distance  at  the 
surface  of  —2900  km.  For  further  technical  details,  see 
section  4.1.1  of  Kidder  and  Vonder  Haar  (1995)  and 
section  3.1  and  appendix  J.l  of  Goodrum  et  al.  (2000). 

This  intrinsic  data  resolution,  known  as  local  area 
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coverage  (LAC),  is  too  dense  to  be  continuously  stored 
and  telemetered  to  ground  stations.  Thus,  LAC  data 
can  only  be  received  at  certain  scheduled  times.  Con¬ 
tinuous  monitoring  at  all  other  times  is  provided  by 
global  area  coverage  (GAC)  data,  which  are  onboard 
averages  of  a  reduced  subset  of  the  raw  LAC  data  as 
described  in  section  10.2.1  of  Kidder  and  Vonder  Haar 
(1995)  and  section  3. 1.3.2  of  Goodrum  et  al.  (2000). 
Surface  footprint  diameters  for  GAC  data  are  ~4  X  3.3 
km^  at  nadir. 

Data  from  each  channel  are  converted  on  board  into 
10-bit  binary  earth  scene  counts  prior  to  being  tele¬ 
metered  to  ground  stations.  Earth  scene  radiances 
are  derived  using  the  nonlinear  radiance  correction  for¬ 
mula, 

~  ^0  ^2^E^  (4) 

where  Uq,  a^,  and  ^2  are  regression  coefficients  derived 
from  prelaunch  calibration  data  (Walton  et  al.  1998). 

b.  Radiosonde  data 

A  meteorological  station  on  Jan  Mayen  launches  ra¬ 
diosondes  at  0000  and  1200  UTC  each  day:  see  Fig.  1  of 
Gjevik  and  Marthinsen  (1978)  for  its  precise  location 
on  the  island.  Here  we  use  winds  and  temperatures 
acquired  from  the  1200  UTC  sounding  on  25  January 
2000  in  our  model  simulations,  which  we  reinterpolate 
onto  a  high  vertical  resolution  grid,  vertically  smooth 
using  a  running  average  of  2-km  width,  then  reinterpo¬ 
late  again  onto  the  models’  vertical  grids.  While  we 
assume  here  that  these  data  approximate  the  true  ver¬ 
tical  profile  over  the  island  at  this  time,  the  actual  bal¬ 
loon  trajectory  is  oblique:  given  a  typical  ascent  velocity 
of  5  m  s“^  (e.g..  Lane  et  al.  2000),  passive  advection 
calculations  based  on  radiosonde  winds  for  this  day 
place  the  balloon  —100  km  downstream  by  the  time  it 
reached  10-km  altitude. 

c.  Meteorological  analyses 

To  provide  some  cross  validation  for  the  model  runs, 
we  also  use  profiles  of  winds  and  temperatures  at  the 
closest  grid  point  to  Jan  Mayen  from  the  12-hourly  Met 
Office  (UKMO)  3.75°  X  2.5°  analyses  (Swinbank  and 
O’Neill  1994).  To  define  surface  weather  patterns,  we 
use  6-hourly  analyzed  mean  sea  level  pressures  (MSLP) 
from  the  National  Centers  for  Environmental  Predic¬ 
tion-National  Center  for  Atmospheric  Research 
(NCEP-NCAR)  2.5°  X  2.5°  reanalyses  (Kalnay  et  al. 
1996). 


4.  Jan  Mayen  wave  clouds  of  25  January  2000 

a.  AVHRR-3  observations 

Figure  3  plots  a  chronological  sequence  of  selected 
AVHRR-3  channel-5  (—11.5-12.5  jixm)  radiances  on  25 
January  2000.  The  data  have  been  reinterpolated  to 
Cartesian  coordinates  {x  being  east-west),  with  the  is¬ 
land  of  Jan  Mayen  centered  at  x  =  y  =  0  on  each  map 
and  its  coastline  plotted  in  each  panel.  We  use  only  the 
thermal  IR  data  since  Jan  Mayen  is  in  polar  night  at  this 
time  of  year. 

At  —1140  UTC,  Fig.  3a  shows  that  the  island  was 
obscured  by  cloud  decks.  Around  3  h  later  at  1500  UTC 
these  clouds  have  been  advected  to  the  east  to  reveal 
evidence  of  mountain-wave  banding  of  lower-level 
clouds:  at  this  time  we  have  both  conventional  GAC 
data  (Fig.  3b)  and  high-resolution  LAC  data  (Fig.  3c). 
The  LAC  image  shows  cloud  banding  superficially  con¬ 
sistent  with  a  purely  diverging  three-dimensional  ship 
mountain-wave  pattern  (Sharman  and  Wurtele  1983) 
emanating  from  Jan  Mayen,  though  only  one  “arm”  of 
this  V-shaped  pattern  is  fully  visible,  the  other  still  be¬ 
ing  partially  obscured  by  overlying  cloud  decks. 

Around  90  min  later  at  1630  UTC,  GAC  radiances  in 
Fig.  3d  again  show  a  diverging  wavelike  cloud  pattern. 
Though  this  is  a  lower-resolution  image  than  the  LAC 
image  in  Fig.  3c,  it  provides  our  least-obscured  image  of 
the  full  wave  pattern.  In  fact  this  wave  pattern  extends 
farther  downstream  in  this  image  than  might  be  sug¬ 
gested  by  the  limited  domain  plotted  in  Fig.  3d. 

Figure  3e  plots  another  AVHRR-3  GAC  image  ac¬ 
quired  roughly  2  h  later  at  —1820  UTC.  Upper-level 
cloud  decks  have  moved  in  and  partially  obscured  the 
wave  pattern.  Nonetheless,  we  still  see  evidence  of  the 
lower  (southward)  arm  of  the  diverging  wave  pattern  in 
the  bottom  half  of  the  image. 

Figure  3f  and  a  sequence  of  later  images  (not  shown) 
do  not  show  evidence  of  banded  clouds  downstream  of 
Jan  Mayen.  While  some,  like  Fig.  3f,  are  heavily  im¬ 
pacted  by  obscuring  cloud  layers,  the  weight  of  evi¬ 
dence  from  these  images  suggests  that  the  banded  wave 
clouds  have  disappeared  at  these  later  times. 

Thus,  from  the  AVHRR-3  channel-5  IR  imagery  in 
Fig.  3,  we  deduce  an  approximate  6-h  duration  of  this 
apparent  Jan  Mayen  wave  event,  lasting  from  —1200  to 
1800  UTC.  The  role  of  high-cloud  decks  in  obscuring 
the  wave  banding  patterns  also  indicates  that  the 
banded  wave  clouds  occur  in  the  lower  or  midtropo¬ 
sphere. 

b.  Synoptic  surface  meteorology 

Figure  4  plots  6-hourly  NCEP-NCAR  reanalyzed 
MSLP  in  a  broad  region  centered  over  Jan  Mayen  from 
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(a)  Jan.  25,  2000:  1 1 40-1 1 43  UTC  (d)  Jan.  25,  2000:  1 639-1 640  UTC 


(b)  Jan.  25,  2000:  1459-1501  UTC  (e)  Jan.  25,  2000:  1819-1820  UTC 
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Fig.  3.  (a)-(f)  Chronological  time  sequence  of  AVHRR-3  channel-5  earth  scene  radiance  images  on  25  Jan  2000  centered  at  x  = 
y  =  0  over  Jan  Mayen,  whose  coastline  is  outlined  in  each  panel.  All  images  are  GAC  data,  except  for  (c),  which  is  LAC  data. 
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Fig.  4.  MSLP  (in  hPa)  derived  from  NCEP-NCAR  reanalysis  fields  at  (a)  0600,  (b)  1200,  (c)  1800  UTC  25  Jan  2000,  and  (d)  0000 
UTC  26  Jan  2000.  Jan  Mayen  is  plotted  as  the  tiny  island  coastline  near  the  center  of  each  map  at  71°N,  8°W. 


0600  UTC  on  25  January  to  0000  UTC  on  26  January. 
We  see  rapid  eastward  passage  of  a  compact  low  pres¬ 
sure  cell  across  the  domain.  As  the  core  of  this  low 
passed  to  the  north  of  Jan  Mayen  at  —1200  UTC,  Fig. 
4b  shows  that  it  produced  enhanced  surface  geostrophic 
westerlies  across  Mount  Beerenberg.  These  surface 
westerlies  persist  over  Jan  Mayen  at  1800  UTC,  though 
the  isobars  in  Fig.  4c  reveal  that  the  core  of  the  low  has 
now  passed  well  to  the  east  of  the  island,  and  that  sur¬ 
face  westerly  flow  across  Mount  Beerenberg  is  about  to 


weaken  and  be  replaced  by  weaker  northerly  flow  as 
the  low  moves  farther  eastward.  This  is  confirmed  by 
the  isobars  at  0000  UTC  in  Fig.  4d. 

Thus,  the  appearance  of  cloud  banding  from  1200 
to  1800  UTC  in  the  AVHRR-3  radiances  in  Fig.  3 
correlates  with  a  period  of  enhanced  surface  westerlies 
across  Mount  Beerenberg  (and  hence  enhanced  po¬ 
tential  for  mountain-wave  forcing)  accompanying  the 
eastward  passage  of  a  polar  low  to  the  north  of  Jan 
Mayen. 
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Fig.  5.  Horizontal  wind  profile  over  Jan  Mayen  at  1200  UTC  25  Jan  2000  from  a  routine 
radiosonde  sounding  from  the  island  (black)  and  from  the  UKMO  analyzed  winds  for  this  day 
and  time  in  its  nearest  3.75°  X  2.5°  grid  box  (gray).  Solid  lines  show  three-dimensional  wind 
profiles,  dashed  lines  show  various  two-dimensional  projections.  Both  profiles  were  smoothed 
by  first  interpolating  the  raw  profile  data  onto  a  higher-resolution  regular  height  grid  and  then 
smoothing  using  a  2-km  sliding  vertical  average. 


c.  Meteorological  profiles  over  Jan  Mayen  at 
1200  UTC 

Figure  5  plots  the  three-dimensional  horizontal  wind 
profile  over  Jan  Mayen  from  the  routine  1200  UTC 
radiosonde  ascent  on  25  January  (solid  curve),  as  well 
as  the  nearest  gridbox  1200  UTC  profile  from  the 
UKMO  analysis.  The  radiosonde  profiles  verify  the 
presence  of  strong  near-westerly  surface-level  winds  in¬ 
ferred  from  Fig.  4b,  with  speeds  —30  m  s“^.  These  wind 
speeds  increase  with  altitude  to  a  jet  maximum  of  —70 
m  s“^  near  10-km  altitude.  This  maximum  is  nearer  55 
m  s“^  in  the  UKMO  profile. 

To  produce  mountain-wave  activity  so  far  down¬ 
stream  of  Jan  Mayen  in  Fig.  3,  it  is  likely  that  these 
waves  were  vertically  trapped  (e.g.,  Simard  and  Peltier 
1982).  Strongly  sheared  wind  profiles  like  those  in  Fig. 
5  are  one  way  to  produce  the  vertical  wave  reflection 
required  for  trapping.  To  assess  this,  we  compute  pro¬ 
files  of  the  Scorer  parameter.  The  standard  Scorer  pa¬ 


rameter  is  a  two-dimensional  term  in  which  the  hori¬ 
zontal  wavenumber  vectors  and  the  wind  vector  are 
coaligned.  For  three-dimensional  problems,  a  range  of 
angles  for  the  horizontal  wavenumber  vector  exist  and 
the  wind  vector  can  rotate  with  altitude,  and  thus  no 
unique  Scorer  parameter  exists.  Some  insights  can  be 
gained,  however,  by  studying  a  generalized  form  for 
three-dimensional  problems. 


Ciz,  <p) 


N\z) 

U\z) 


(5) 


where  N{z)  is  the  background  buoyancy  frequency  pro¬ 
file,  U{z)  =  Utot(z)  cos[q:(z)  -  (p]  is  the  profile  of  the 
component  of  the  wind  vector  f/tot(cos  a,  sin  a)  pro¬ 
jected  along  a  given  horizontal  wavenumber  vector’s 
direction  cp  =  arctan  (//k),  is  the  curvature  in  this 
projected  component  wind  profile,  and  = 

[U^{z)  +  V^{z)Y'^.  An  equivalent  expression  is  utilized 
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Profiles  over  Jan  Mayen:  25  January  2000  1200  UTC 


(a)  Horizontal  Wind:  a=20 


(b)  Horizontal  Wind  Curvature:  a=20° 


m  s  km 

(d)  Scorer  Parameters 


Fig.  6.  Profiles  of  (a)  horizontal  wind  and  (b)  wind  curvature  projected  along  an  axis  directed  a  =  20°  north  of  east,  (c)  Brunt-Vaisala 
(buoyancy)  frequencies,  (d)  Scorer  parameter  profiles  cp)  for  cp  =  a  =  20°  with  curvature  omitted  (thick  solid)  and  curvature 

retained  (dotted  curves).  Dashed  curves  show  €^(z,  9)  -  with  curvature  omitted  for  A/^  =  277//:;,  =  30  km.  In  all  cases,  black  curves 
are  from  the  Jan  Mayen  radiosonde  ascent  and  gray  curves  are  from  UKMO  analysis  at  1200  UTC  25  Jan  2000.  All  profiles  are 
smoothed  vertically  by  resampling  at  high  vertical  resolution  and  then  smoothing  using  a  2-km  running  average,  except  for  the  thin 
curve  in  (a),  which  plots  the  raw  radiosonde  winds. 


by  Vosper  and  Worthington  (2002)  and  Eckermann  et 
al.  (2006).  From  the  profiles  in  Fig.  5,  specifically  the 
hodographic  projection  at  the  surface,  we  see  that  the 
horizontal  wind  vector  is  directed  at  an  approximately 
constant  angle  of  a  =  20°  north  of  due  east  through 
most  of  the  troposphere. 

Figures  6a, b  plot  component  profiles  of  wind  speed 
and  curvature,  respectively,  for  both  the  radiosonde 
and  UKMO  profiles  on  25  January  at  1200  UTC,  using 
a  choice  for  the  direction  of  the  horizontal  wavenumber 
of  (p  =  a  =  20°  (i.e.,  wave  vectors  aligned  parallel  to  the 
wind  vector).  Figure  6c  plots  the  buoyancy  frequency 
profiles.  The  thin  black  curve  in  Fig.  6a  shows  the  raw 


radiosonde  wind  profile,  while  the  thick  black  curve 
shows  that  same  profile  after  high-resolution  resam¬ 
pling  and  then  smoothing  using  a  sliding  2-km  vertical 
averaging  window  to  remove  the  small-scale  structure. 
This  filter  is  applied  to  all  the  profiles  (both  radiosonde 
and  UKMO)  that  enter  into  the  Scorer  parameter  cal¬ 
culation  (5),  since  even  slight  vertical  structure  in  the 
wind  and  temperature  profiles  gets  greatly  enhanced 
during  the  numerical  gradient  and  curvature  calcula¬ 
tions  in  (5),  yielding  large  spurious  oscillations  in  £^(z, 
(p)  (see,  e.g.,  Danielsen  and  Bleck  1970;  Shutts  1992; 
Smith  2004). 

The  thick  curves  in  Fig.  6d  show  Scorer  parameters 
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20°)  calculated  by  ignoring  the  curvature  term  in 
(5),  since  the  first  term  N^{z)lu{z)  is  generally  the 
most  significant.  Furthermore,  curvature  terms  are 
omitted  at  present  from  our  Fourier-ray  model’s  wave 
equations,  so  this  simplified  term  is  roughly  equivalent 
to  this  model’s  Scorer  parameter.  (The  Scorer  param¬ 
eter  and  the  present  version  of  the  Fourier-ray  model 
omit  Coriolis  terms,  but  their  influence  on  these  short- 
wavelength  trapped  wave  solutions  is  minimal.)  The 
corresponding  dotted  curves  show  20°)  with  the 
curvature  profiles  in  Fig.  6b  retained  in  the  calculations. 
They  produce  fairly  minimal  modifications  to  the  sim¬ 
plified  N^(z)IU^(z)  profile.  Tests  that  include  addi¬ 
tional  wind  shear  and  density-scale  height  terms  in  a 
generalized  Scorer  parameter  derived  from  fully  com¬ 
pressible  wave  equations  (Nance  1997)  also  yielded 
minimal  changes  to  these  simplified  N'^{z)lu{z)  pro¬ 
files. 

Figure  6d  shows  that  the  wind  and  stability  profiles 
over  Jan  Mayen  on  25  January  at  1200  UTC  yield  a 
deep  minimum  in  €^(z,  9)  at  altitudes  —5-10  km.  Turn¬ 
ing  points  Zt  for  a  given  mountain  wave  of  total  hori¬ 
zontal  wavelength  A;,  =  27Tlk^  occur  where  €^(z,  9)  - 
k\  =  0.  The  dashed  curves  in  Fig.  6d  plot  €^(z,  20°)  -  k\ 
for  =  30  km.  We  see  that  the  radiosonde-derived 
Scorer  parameter  predicts  reflection  of  this  horizontal 
wavelength  at  ^  5  km,  whereas  the  UKMO  profile 
predicts  no  reflection  (due  to  its  weaker  peak  winds), 
though  it  comes  very  close  to  yielding  reflection  at  —8 
km.  Since  waves  of  comparable  or  shorter  horizontal 
wavelengths  are  likely  to  be  generated  by  flow  over  the 
horizontally  narrow  obstacle  presented  by  Mount  Beer- 
enberg  (see  Fig.  1),  we  see  that  Jan  Mayen  mountain 
waves  probably  reflected  vertically  in  the  midtropo¬ 
sphere  at  1200  UTC.  Since  we  have  a  stable  boundary 
layer  (Fig.  6c)  containing  fairly  strong  surface  flow  (Fig. 
6a),  these  reflected  waves  will  reflect  again  at  the  sur¬ 
face,  rather  than  being  absorbed  (Smith  et  al.  2002), 
yielding  a  vertically  trapped  wave  that  can  penetrate 
long  distances  downstream  within  this  waveguide,  as 
observed  in  Fig.  3. 

The  surface  wind  speed  Uq  is  ^30  m  s“^  (Fig.  6a)  and 
surface  buoyancy  frequency  Nq  is  ^0.015  rad  s“^  (Fig. 
6c).  Since  the  maximum  height  of  the  Jan  Mayen  to¬ 
pography  is  ^2  km  (Fig.  2),  this  yields  Fr“^  ^  1, 
where  Fr  =  U^Noh^  is  the  surface  Froude  number. 
This  suggests  a  near-maximum  amplitude  for  the  forced 
wave,  with  the  surface  flow  passing  over  (rather  than 
around)  the  mountain  and  translating  all  of  the  moun¬ 
tain  pressure  drag  into  a  gravity  wave  response.  Since 
Fr“^  is  close  to  unity,  some  weak  gravity  wave  breaking 
and/or  associated  lee-vortex  nonlinearity  may  occur  in 
a  thin  layer  near  the  surface. 


5.  Model  results 

a.  Fourier-ray  solutions 

We  begin  by  studying  how  horizontal  cross  sections 
of  our  Fourier-ray  model  solutions  compare  with  the 
AVHRR-3  imagery  in  Fig.  3.  Our  best  (least  obscured) 
AVHRR-3  image  of  the  wavelike  cloud  banding  in  Fig. 
3d  is  replotted  in  Fig.  7a  as  our  observational  reference. 
Below  it,  we  plot  two  Fourier-ray  solutions  for  vertical 
velocity  at  z  =  3  km  and  r  =  4  h,  as  calculated  using  the 
Jan  Mayen  radiosonde  profiles  (Fig.  7b)  and  the 
UKMO  profiles  (Fig.  7c). 

Both  simulations  show  similarities  with  the 
AVHRR-3  image  in  Fig.  7a.  The  major  response  in 
each  case  is  a  V-shaped  trapped  mountain-wave  pattern 
extending  long  distances  downstream  at  angles  to  the 
island  similar  to  those  observed  in  Fig.  7a.  There  are, 
however,  some  noticeable  differences  between  the  two 
solutions. 

Sharman  and  Wurtele  (1983)  showed  that  analytical 
solutions  for  trapped  mountain  waves  from  three- 
dimensional  symmetric  obstacles  have  properties  and 
sensitivities  that  resemble  classical  Kelvin  solutions  for 
ocean  surface  ship  wakes  (Reed  and  Milgram  2002). 
Like  ship  wakes,  trapped  mountain  waves  can  exhibit 
both  transverse  waves,  which  appear  downstream  of 
the  obstacle  as  linear  phase  lines  orthogonal  to  the  flow 
direction,  and  diverging  waves,  which  appear  as  curved 
phase  lines  in  two  V-shaped  wedges  of  activity  either 
side  of  the  downstream  flow  axis. 

The  Fourier-ray  model  simulation  based  on  the 
UKMO  profiles  in  Fig.  7c  contains  both  diverging  and 
transverse  waves,  but  the  solution  based  on  the  Jan 
Mayen  radiosonde  profile  in  Fig.  7b  contains  only  di¬ 
verging  waves.  The  AVHRR-3  images  in  Figs.  3  and  7a 
show  clear  evidence  of  diverging  waves  but  no  evidence 
of  transverse  waves  in  the  cloud  banding:  see  especially 
the  high-resolution  LAC  image  in  Fig.  3c.  This  indicates 
that  the  (presumably)  more  accurate  representation  of 
the  meteorological  profiles  over  Jan  Mayen  from  the 
radiosonde  ascent  is  important  for  an  accurate  simula¬ 
tion  of  the  observed  disturbance.  We  will  have  more  to 
say  in  section  6b  about  why  the  transverse  waves  are 
present  for  the  UKMO  profile  simulations  but  not  for 
the  radiosonde  profile  simulations. 

Other  aspects  of  the  radiosonde-based  Fourier-ray 
solution  also  agree  better  with  the  imagery  than  does 
the  UKMO-based  solution.  For  example,  the  radio- 
sonde-based  solution  shows  greater  “flaring”  of  the  di¬ 
verging  waves  into  progressively  wider  wedge  angles 
with  increasing  distance  downstream.  By  wedge  angle, 
we  mean  the  angle  relative  to  a  line  passing  from  the 
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(c)  Fourier-ray  with  UKMO  profiles 
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Fig.  7.  (a)  Wave  train  detected  in  AVHRR-3  image  on  25  Jan  2000,  enlarged  from 
Fig.  3d.  (b)  The  Fourier-ray  solution  for  vertical  velocity  at  z  =  3  km  and  t  =  4  h  using 
the  radiosonde  wind  and  stability  profiles  in  Fig.  6.  Minimum  value  is  -6.5  m  and 
maximum  value  is  5.9  m  s“^.  (c)  Same  as  in  (b),  but  using  the  UKMO  wind  and  stability 
profiles  in  Fig.  6.  Minimum  value  is  -7.6  m  and  maximum  value  is  7.9  m  s“^.  The 
bars  in  (b)  and  (c)  show  vertical  velocity  values  in  m  s“^. 
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Fig.  8.  Vertical  velocity  at  z  =  3  km  and  t  =  4  ti,  computed  from  the  nonlinear  numerical 
model  using  the  radiosonde  profiles  in  Fig.  6.  The  bar  shows  the  range  in  m  s“^.  Minimum 
value  is  -7.3  m  s“^  and  maximum  value  is  4.2  m  s“^.  The  coordinate  axes  and  the  amplitude 
scale  are  the  same  as  for  the  corresponding  Fourier-ray  solution  in  Fig.  7b.  The  line  that  slopes 
upward  from  left  to  right  indicates  the  horizontal  locus  of  the  vertical  cross  section  in  Fig.  9b. 
The  line  makes  an  angle  of  36°  with  the  x  axis. 


center  of  the  topography  through  the  approximate  cen¬ 
ter  of  the  downstream  wave  field.  There  is  also  a  clear 
horizontal  wavelength  dependence  to  the  diverging 
waves,  with  the  shortest  wavelengths  confined  to  the 
smaller  wedge  angles  and  the  longer  wavelengths  ex¬ 
tending  into  larger  wedge  angles.  A  similar  wavelength 
dependence  is  seen  in  idealized  numerical  model  solu¬ 
tions  (see,  e.g..  Fig.  9a  of  Sharman  and  Wurtele  2004). 

Since  the  solutions  based  on  the  radiosonde  profile 
agree  better  with  the  AVHRR-3  observations,  we  focus 
on  them  for  further  analysis. 

b.  Nonlinear  numerical  model  simulation 

To  assess  more  rigorously  the  accuracy  of  the  Fou¬ 
rier-ray  solutions  in  Fig.  7,  we  performed  companion 
runs  using  a  nonlinear  numerical  model  containing 
identical  topography,  as  described  in  sections  2b  and  2c. 
Figure  8  plots  this  model’s  vertical  velocity  field  at  z  = 
3  km  after  4  h  using  the  radiosonde  profiles.  The  line 
sloping  upward  from  left  to  right  indicates  the  horizon¬ 
tal  coordinates  of  a  vertical  cross  section  that  is  shown 
in  Fig.  9b. 

The  similarities  between  the  numerical  model  solu¬ 
tion  in  Fig.  8  and  the  corresponding  Fourier-ray  solu¬ 
tion  in  Fig.  7b  are  striking.  Both  solutions  reproduce  a 
purely  diverging  trapped  wave  solution  with  similar 


wavelengths  and  phase  orientation.  Both  also  include 
very  short  wavelengths  on  the  inside  of  the  upper  arm 
of  the  V-shaped  pattern,  located  just  below  the  line 
indicating  the  position  of  the  vertical  cross  section,  with 
the  phases  nearly  parallel  to  that  line. 

Figure  9  plots  height-varying  cross  sections  of  the 
vertical  velocities  along  the  line  in  Fig.  8  for  the  Fou¬ 
rier-ray  model  (Fig.  9a)  and  the  nonlinear  numerical 
model  (Fig.  9b).  Both  show  a  disturbance  trapped  be¬ 
tween  the  ground  and  ~8  km  altitude  with  similar  hori¬ 
zontal  wavelength  structure  between  0  and  8  km.  Dif¬ 
ferences  at  higher  altitudes  and  downwind  are  dis¬ 
cussed  below. 

6.  Discussion 

Our  Fourier-ray  simulations  of  mountain  waves 
forced  by  observed  flow  over  Jan  Mayen  on  25  January 
2000  produced  vertical  velocity  fields  iv(x,  y,  z,  r)  at  z  = 
3  km  and  ^  =  4  h  that  reproduced  the  salient  features  of 
high-resolution  banded  IR  cloud  imagery  acquired 
from  space  by  AVHRR-3.  For  further  validation,  we 
also  compared  those  model  results  with  output  from 
companion  runs  using  a  nonlinear  numerical  model, 
again  revealing  good  overall  agreement. 

On  a  2.8-GHz  single  processor  workstation,  these 
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(a)  Fourier-Ray  Solution 
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Fig.  9.  Vertical  cross  sections  of  vertical  velocity  in  m  s“^  (see 
grayscale  bars)  at  t  =  4  h  along  the  line  indicated  in  Fig.  8.  (a) 
Fourier-ray  solution.  Minimum  value  is  -7.8  m  s“^  and  maximum 
value  is  6.1  m  s“^.  (b)  Nonlinear  numerical  model  solution.  Mini¬ 
mum  value  is  -8.7  m  and  maximum  value  is  5.8  m  s“^. 


Fourier-ray  solutions  took  ~5  min  to  compute  in 
FORTRAN  90  on  a  1024  X  1024  horizontal  grid  at  40 
levels  from  the  ground  to  20  km.  The  corresponding 
nonlinear  numerical  model  runs  at  256  X  256  X  32  with 
a  4-s  time  step  on  the  same  workstation  took  38  h  to 
yield  the  ^  =  4  h  solution.  Since  the  Lipps-Hemler  code 
we  use  has  not  been  aggressively  optimized,  we  per¬ 
formed  identical  runs  on  the  same  workstation  using  the 
Weather  Research  and  Forecasting  (WRF)  Advanced 
Research  (ARW)  model  (Skamarock  et  al.  2005). 
These  ARW  runs  took  21  h  to  yield  the  ^  =  4  h  solution 
using  a  4-s  time  step,  and  8.5  h  using  a  10-s  time  step. 


Some  differences  among  the  various  model  solutions 
were  noted.  We  investigate  them  further  here. 

a.  Vertical  wave  structure 

The  vertical  cross  sections  of  the  Fourier-ray  and  nu¬ 
merical  model  solutions  in  Fig.  9,  while  similar  in  some 
respects,  also  show  some  differences.  Trapped  waves  in 
the  nonlinear  numerical  model  solution  decay  in  am¬ 
plitude  more  rapidly  with  downstream  distance  than 
the  Fourier-ray  solution,  while  at  stratospheric  altitudes 
the  nonlinear  numerical  model  solution  contains  more 
energy  in  freely  propagating  (phase  tilted)  mountain 
waves  than  the  Fourier-ray  solution. 

Both  differences  seem  to  be  related  to  the  tunnelling 
of  trapped  waves  (Sutherland  and  Yewchuk  2004) 
through  an  evanescent  layer  surrounding  the  wind  jet  at 
heights  —7-10  km,  which  then  emerge  into  a  free- 
propagating  region  above  —10  km  (see  Fig.  6d).  This 
“leakiness”  of  the  tropospheric  duct  leads  to  greater 
decay  with  downwind  distance  of  the  vertically  trapped 
tropospheric  wave  amplitudes  in  Fig.  9b.  Similar  tun¬ 
nelling  of  trapped  mountain  waves  into  the  strato¬ 
sphere  was  noted  by  Eckermann  et  al.  (2006)  in  a 
Lipps-Hemler  model  forecast  for  this  same  date  over 
northern  Scandinavia,  a  region  —30°  to  the  east  of  this 
Jan  Mayen  wave  event. 

Our  Fourier-ray  model  does  not  yet  account  for  tun¬ 
nelling,  which  requires  a  two  turning-point  analysis 
(Bender  and  Orszag  1978).  One  turning  point  is  at  the 
lower  edge  of  the  evanescent  region,  and  the  other  is  at 
the  upper  edge.  Our  single  turning-point  Fourier-ray 
solution  breaks  down  near  and  above  the  upper  turning 
point.  Here,  as  a  temporary  measure,  we  have  set  the 
wave  amplitude  to  zero  whenever  the  waves  approach 
the  second  turning  point.  Specifically,  we  set  the  wave 
amplitude  to  zero  whenever  z  >  Ztp  +  O-^Az^p,  where 
Ztp  is  the  height  of  the  lower  turning  point  and  Aztp  is 
the  height  difference  between  the  upper  and  lower 
turning  points.  This  was  done  only  for  waves  that  have 
two  turning  points.  The  solutions  for  the  vertically 
propagating  waves  and  for  waves  with  a  single  turning 
point  were  left  unchanged. 

To  include  the  effects  of  tunnelling  in  the  Fourier-ray 
model,  we  would  need  either  to  match  the  ray  solution 
to  separate  Airy  function  solutions  around  each  of  the 
two  turning  points,  or  to  introduce  a  uniformly  valid 
two  turning-point  solution  in  terms  of  parabolic  cylin¬ 
der  functions  [sometimes  called  Weber  functions:  see 
Kravtsov  and  Orlov  (1999)].  We  hope  to  test  such  so¬ 
lutions  in  future  work.  We  have,  however,  estimated 
how  much  attenuation  should  occur  through  the  eva¬ 
nescent  layer  in  the  current  problem  by  evaluating 
exp(-/^  \m\  dz)  for  each  k,  I  value  that  results  in  two 
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turning  points.  Here  a,  b  are  the  heights  of  the  lower 
and  upper  turning  points,  respectively,  and  the  expo¬ 
nential  is  a  measure  of  the  wave  attenuation  in  the 
evanescent  region  between  the  turning  points.  We 
found  that  important  parts  of  the  spectrum  had  attenu¬ 
ation  factors  greater  than  0.1,  suggesting  that  their  tun¬ 
nelled  wave  amplitudes  would  be  noticeable  above  the 
wind  jet  maximum.  For  example,  for  the  k,  I  value  con¬ 
sidered  in  calculating  the  Scorer  parameter  of  Fig.  6d, 
the  value  of  the  exponential  is  about  0.6. 

b.  Resonant  modes 

While  the  radiosonde  profile  produced  a  purely  di¬ 
verging  trapped  wave  solution  in  Fig.  7b,  the  UKMO 
profile  produced  a  solution  with  both  diverging  and 
transverse  components  (Fig.  7c). 

To  investigate  this  further,  we  calculated  the  reso¬ 
nant  modes  for  each  of  these  solutions.  These  are  the 
singularities  of  the  steady-state  {t  oo)  trapped  wave 
solution  (All),  which  involve  the  Airy  function  Ai(r). 
The  Airy-function  argument  r,  defined  in  (A8),  is  posi¬ 
tive  in  the  evanescent  region  above  the  turning  point, 
negative  in  the  propagation  region  below  the  turning 
point,  and  zero  at  the  turning  point.  Resonances  occur 
whenever  Ai(ro)  =  0,  with  Vq  =  r{z  =  0).  This  occurs 
approximately  at  Vq  =  —2.34,  —4.09,  and  —5.52,  for 
modes  1,  2,  and  3,  respectively. 

Contours  of  Vq  as  a  function  of  the  horizontal  wave- 
numbers  k,  I  are  plotted  in  Fig.  10  for  the  radiosonde 
profiles  (Fig.  10a)  and  for  the  UMKO  profiles  (Fig. 
10b).  The  bold  curves  indicate  the  resonant  values 
where  Ai(ro)  =  0.  A  gap  in,  or  termination  of,  a  par¬ 
ticular  Eq  contour  occurs  where  the  local  k,  I  values  do 
not  correspond  to  a  trapped  wave.  Generally,  the  k,  I 
values  to  the  upper  right  of  the  plotted  contours  corre¬ 
spond  to  propagating  waves  without  turning  points, 
while  the  k,  I  values  to  the  lower  left  of  the  plotted 
contours  correspond  to  evanescent  waves  whose  intrin¬ 
sic  frequency  is  greater  than  N  at  the  ground. 

For  the  radiosonde  profiles,  there  is  a  gap  in  the 
mode-1  resonance  curve,  centered  around  ktk^  =  -0.5. 
For  the  UKMO  case,  the  mode-1  curve  is  continuous. 
This  seems  to  be  the  reason  why  transverse  waves  occur 
for  the  UKMO  profile  solution  in  Fig.  7c  but  not  for  the 
radiosonde  profile  solution  in  Fig.  7b. 

Spatial-ray  tracing  illustrates  the  point.  Figure  11 
plots  horizontal  group  trajectories  of  some  rays  with  {k, 
1)  values  sampled  from  the  corresponding  mode-1  reso¬ 
nance  curves  in  Fig.  10.  These  trajectories  were  com¬ 
puted  by  numerically  integrating  the  spatial  ray  equa¬ 
tions  (Lighthill  1978) 

(6) 


(a)  Radiosonde  Profile 


Fig.  10.  Contours  of  r^,  the  argument  of  the  Airy  function  so¬ 
lution  in  (All),  for  (a)  the  radiosonde  profile  and  (b)  the  UKMO 
profile.  The  bold  lines  are  the  resonant  values  and  are  labeled  by 
their  mode  number,  1, 2,  or  3,  corresponding  to  Tq  =  —2.34,  —4.09, 
and  -5.52,  respectively.  The  other  Tq  contours  are  separated  at 
intervals  of  -1,  starting  with  Tq  =  -1  (lowermost-left  contour). 
The  horizontal  wavenumbers  k,  I  on  the  axes  are  normalized  (ar¬ 
bitrarily)  by  =  lirlilt)  km). 

dm/dt  =  -kU^  -  IV^  -  coN^/N,  (7) 

out  to  4  h,  the  same  time  at  which  the  corresponding 
Fourier-ray  solutions  in  Figs.  7b,c  were  evaluated.  Here 
X  is  the  position  vector  of  the  ray  and  is  its  ground- 
based  group  velocity  vector. 

Figure  11  shows  that  the  presence  (UKMO  case)  or 
absence  (radiosonde  case)  of  mode-1  rays  in  the  region 
directly  downwind  of  the  mountain  corresponds  to  the 
presence  or  absence  of  transverse  waves  in  the  Fourier- 
ray  solutions  in  Fig.  7.  For  the  radiosonde  profile,  there 


dx/dt  =  Cg, 
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(a)  Radiosonde  Profile 


(b)  UKMO  Profile 


Fig.  11.  The  horizontal  coordinates  of  ray  paths  corresponding 
to  points  on  the  mode-1  resonance  curves  in  Fig.  10.  (a)  For  the 
radiosonde  profile,  66  ray  paths  are  shown,  with  klk^  values  in  the 
range  of  -1.5  to  -0.25.  (b)  For  the  UKMO  profile,  50  ray  paths 
are  shown,  with  /cZ/cq  values  in  the  range  of  —1.4  to  —0.88. 

are  rays  with  wavenumbers  {k,  1)  lying  off  of  this 
mode-1  curve  that  propagate  downwind.  However, 
these  rays  are  either  vertically  propagating  (as  opposed 
to  vertically  trapped),  and  thus  unable  to  produce  a 
transverse  wave  train  at  long  downwind  distances,  or 
they  are  vertically  trapped  but  nonresonant.  Resonance 
is  important  because  these  are  the  waves  that  construc¬ 
tively  interfere  in  the  solution  (A3)  to  produce  the  larg¬ 
est  downstream  amplitudes.  We  also  checked  the  rays 
for  the  mode-2  resonances  in  the  radiosonde  data,  and 
again  we  found  an  absence  of  rays  in  this  same  region. 

Note  that  the  spatial-ray  tracing  also  provides  other 
consistency  checks  on  our  Fourier-ray  solutions.  For 
example.  Fig.  11  shows  rays  governed  by  the  radio¬ 
sonde  profile  “flaring”  more  in  the  diverging  wings  and 
propagating  farther  downstream  after  4  h  than  those 


rays  governed  by  the  UKMO  profile,  consistent  with 
the  Fourier-ray  solutions  in  Fig.  7. 

The  physical  interpretation  of  these  results  is  as  fol¬ 
lows.  Transverse  waves  are  absent  in  the  radiosonde- 
based  solution  because  its  mode-1  resonant  wavenum¬ 
ber  magnitudes  in  Fig.  10  are  smaller  than  those  of 
the  UKMO-based  solution.  The  smaller  radiosonde- 
based  mode-1  kh  values  do  not  satisfy  the  reflection 
criterion  €^(z,  cp)  -  =  0  for  the  transverse  wavenum¬ 

ber  alignments  9,  and  so  these  rays  propagate  freely  in 
the  vertical,  leaving  a  purely  diverging  trapped  wave 
pattern.  Conversely,  the  larger  k^  values  on  the 
UKMO-based  resonant  mode-1  curve  satisfy  €^(z,  9)  - 
=  0  for  all  wavenumber  alignments,  resulting  in  a 
trapped  wave  solution  with  both  diverging  and  trans¬ 
verse  components.  The  smaller  kf^  (longer  resonant 
values  of  the  radiosonde-based  solution  compared  to 
the  UKMO-based  solution  are  evident  in  the  wave  pat¬ 
terns  in  Figs.  7b,c. 

7.  Summary 

We  have  used  the  Fourier-ray  method  to  simulate 
mountain  waves  that  were  observed  by  high-resolution 
IR  satellite  cloud  imagery  and  that  extended  long  dis¬ 
tances  downwind  of  the  island  of  Jan  Mayen  on  25 
January  2000.  This  was  a  first  step  in  applying  the  Fou¬ 
rier-ray  method  to  a  real-world  problem.  The  results 
highlight  some  of  the  desirable  features  of  the  Fourier- 
ray  method,  such  as  the  ability  to  produce  very  high 
resolution  solutions  quickly  for  realistic  topography 
and  arbitrary  wind  and  stability  profiles. 

Comparison  with  a  companion  simulation  using  a 
nonlinear  numerical  model  of  flow  over  orography  re¬ 
vealed  similar  wave  patterns  but  also  indicated  the  pos¬ 
sible  importance  (see  Fig.  9)  of  wave  tunnelling  through 
the  wind  jet  centered  at  a  height  of  about  10  km.  Wave 
tunnelling  (i.e.,  the  penetration  of  waves  through  a  lo¬ 
calized  evanescent  region)  is  a  two  turning-point  effect, 
which  is  not  included  in  the  present  formulation  of  the 
Fourier-ray  model.  The  analysis  of  two  turning-point 
problems  is  standard  (Bender  and  Orszag  1978),  and 
we  hope  to  add  this  capability  to  the  Fourier-ray  model 
in  the  future. 

The  speed  and  accuracy  of  this  algorithm,  coupled 
with  its  use  of  arbitrary  surface  topography  and  vertical 
meteorological  profiles,  make  it  promising  as  a  moun¬ 
tain-wave  forecasting  tool,  particularly  as  a  potential 
next-generation  dynamical  core  for  the  MWFM.  Future 
work  will  apply  the  method  to  more  complex  real-world 
problems,  such  as  extended  topography  and  profiles 
based  on  operational  NWP  fields  (see  Eckermann  et  al. 
2004). 
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APPENDIX 

The  Transient  Fourier-Ray  Solution 


where  p  is  the  background  density  and  c^3  is  the  vertical 
group  velocity.  Here  Gq  is  G  at  z  =  0. 

The  wave  phase  is 

(/)=-  m{kj,z')dz'.  (A6) 

J  0 

The  sign  convention  here  is  m  >  0,  with  a  negative  sign 
inserted  in  (A6)  to  produce  downward  phase  propaga¬ 
tion  and  upward  group  propagation. 

For  the  vertically  trapped  waves,  Ai  is  the  Airy  func¬ 
tion,  and 


We  consider  mountain  waves  that  begin  to  emerge 
from  the  mountain  at  ^  =  0  and  are  continuously  gen¬ 
erated  thereafter.  We  assume  that  all  rays  have  the  full 
wave  amplitude  of  the  longtime  steady-state  solution. 
The  resulting  time-dependent  ray  solutions  were  used 
to  produce  the  results  in  the  present  study  and  are  given 
below.  They  were  derived  in  Broutman  et  al.  (2006), 
who  expressed  them  in  terms  of  the  vertical  displace¬ 
ment  r)(k,  I,  z,  t)-  This  is  related  to  the  vertical  velocity 
(used  in  the  present  study)  by  w  =  — /cof),  where  o)  is  the 
intrinsic  frequency. 

The  ray  solution  is  the  sum  of  solutions  for  the  ver¬ 
tically  propagating  waves  and  the  vertically  trapped 
waves: 


^  =  1(4>^  +  4>2),  (A7) 


r  = 


(A8) 


The  phases  and  4>2  are  for  a  ray  incident  upon  and 
reflected  from  the  caustic,  respectively. 

The  time  dependence  of  the  trapped  wave  solution 
(A3)  comes  in  through  defined  by 

^  — /MO 


where 


=  Tjpr  +  'Htr-  (Al) 

For  the  vertically  propagating  waves 

fjp,  =  Fh[G^Gf'^e*t  (A2) 

For  the  vertically  trapped  waves 

Tj,,  =  2nT^'^h[Go/Gf'^[-rf^Ai{r)e‘^^-^'^^Sj^t.  (A3) 

These  are  Eqs.  (5)  and  (22),  respectively,  of  Broutman 
et  al.  (2006),  with  the  latter  equation  assuming  equal 
numbers  of  upgoing  and  downgoing  waves.  The  terms 
in  these  equations  are  defined  below. 

The  factor  F  accounts  for  the  time  dependence  in  the 
vertically  propagating  waves  and  is  defined  by 

Cg3(k,l,z,t')dt'  -Z^j.  (A4) 

The  Heaviside  function  N(^)  is  zero  for  ^  <  0  and  unity 
for  g  >  0.  The  integral  in  the  argument  of  H  is  the 
height  at  time  t  of  the  initial  ray  (for  that  k,  1)  generated 
at  the  mountain  at  ^  =  0. 

The  topography  h{x,  y)  has  a  Fourier  transform  h{k,  /). 
The  factor  G,  defined  such  that  G|fiP  is  the  vertical 
flux  of  wave  action,  is 

G  =  pN^Cg^/d),  (A5) 


O  =  -2 


m{k,  I,  z')  dz'  —  77/2. 


(AlO) 


The  counting  variable  M  indicates  the  number  of  pairs 
of  incident  and  reflected  rays.  Initially  M  =  0  with  5*0 
defined  to  be  zero.  Each  time  the  initial  ray  makes  a 
complete  round  trip  from  the  ground  to  the  turning 
point  and  back  to  the  ground,  M  increases  by  1.  We 
compute  M  as  a  function  of  time  using  Eq.  (33)  of 
Broutman  et  al.  (2006). 

For  the  vertically  trapped  waves,  the  longtime 
steady-state  solution  is 


fjtj-(/c,  l,z,t^^)  =  h 


Gq 

1/2 

r 

_G  _ 

Jo. 

Ai(r) 
Ai(ro)  ’ 


(All) 

where  Eq  is  r  at  z  =  0.  This  is  Eq.  (A9)  of  Broutman  et 
al.  (2003).  It  is  used  here  in  section  6b  for  the  calcula¬ 
tion  of  resonant  modes,  which  occur  where  Ai(ro)  =  0. 
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